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Abstract 

The positronium molecule (PS2) has not been experimentally observed yet 

because its tiny (4.5 eV) binding energy cannot be detected when the molecule 

annihilates by emitting two photons with energy of 0.51 MeV each. It is shown 

in this paper that the electric dipole transition between the recently found 

L = 1 excited-state and the L = ground-state with its characteristic photon 

energy of 4.94 eV is a clear signature of the existence of the positronium 

molecule and the possibility of its experimental observation is realistic. The 

probability of this transition is about 17 % of the total decay rate. An other 

Coulomb four-body system containing positron, HPs (the positronium hydride 

or hydrogen positride), is also included for comparison. 
PACS numbers: 36.10.Dr,31.15.Pf 
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I. INTRODUCTION 



Despite the early theoretical prediction of its existence M, the PS2 molecule has not 
been experimentally found to date. The difficulty stems from the fact that this system 
is neutral and therefore it cannot be separated from the positronium atoms (Ps) and its 
primary decay mode, the annihilation by two-photon emission, is exactly the same as that 
of the Ps atom. The energy of the photons arising from the annihilation is different in 
principle: The photons carry 1.02 MeV energy due to the annihilation plus the binding 
energy of the corresponding system. The binding-energy difference is, however, less than 
1 eV and adding it to 1.02 MeV, the energy of the photons coming from the Ps atom or 
Ps2 molecule cannot be experimentally distinguished. The experimental observation of the 
biexcitons can be considered as an indirect indication of the existence of Ps2- 

In our recent Letter || we have predicted the existence of a hitherto unknown bound 
excited-state of the Ps 2 molecule. In this paper we give a detailed description of this state. 
We have investigated possible decay modes of this state with a special emphasis on the 
electric dipole (El) transition to the ground state. It will be shown that the probability 
of the El transition is comparable to that of the annihilation. The unique energy of this 
transition may possibly be utilized as a sign for the experimental identification of the Ps2 
molecule. 

The stochastic variational method HQ has been used to solve the Coulomb four-body 
problem. In this method the variational trial functions are optimized by gambling: Ran- 
domly chosen configurations are probed and most adequate functions are selected to be the 
basis states. 

The Correlated Gaussians (CG) || are used as basis functions in this procedure. The CG 
basis has a long history in atomic and molecular physics and highly accurate calculations 
are based on this form of basis functions P,|^-[T0[| . The angular part is given by the global 
vector representation |7j . This approach greatly simplifies the calculations for non-spherical 
systems by replacing the partial wave expansion with a much simpler representation of the 
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angular motion. 

The hydrogen positride (positronium hydride), HPs, has already been in the focus of 
intensive theoretical and experimental investigation. This is an ideal system to test the 
SVM. We compare the properties of the Ps 2 and HPs molecules. 

The plan of this paper is as follows. In sect. II we give a brief description of the trial 
function and the stochastic variational method. In sect. Ill the results are presented. The 
main results of the paper are summarized in sect. IV. In Appendices A-D we collect some 
basic ingredients which are used in the present study in order to help readers reproduce our 
results: formulae of the matrix elements in the CG basis, the separation of the center-of- 
mass motion from the CG basis, the use of Sherman- Morrison formula in selecting nonlinear 
parameters, and the symmetry requirement for the trial wave function of the Ps 2 molecule. 



A system of two electrons with mass m and two positive unit charges of mass M is 
considered. Their relative mass is characterized by the ratio a = m/M, and the positronic 
limit is realized by a = 1. (Though we consider the case of a — 1 in this paper, the extension 
to other a values is straightforward, so we give a formulation assuming an arbitrary mass 
ratio.) The Hamiltonian of the system reads as 



where qi and are the charges and the position vectors of the particles. Particle labels 1 and 
3 denote the positive charges, while labels 2 and 4 denote the negative charges. A relative 
coordinate system is introduced by defining xi and x 2 as the distance vectors between the 
positive and negative charges in the first and second atom, and x 3 as the distance vector 
between the center-of-masses of the two atoms: 



II. THE CALCULATION 



4 

H = J2T t -T cm + J2 

i=l i<j 




(1) 



X! = n - r 2 , 



(2) 



x 2 = r 3 - r 4 , 



(3) 
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Mr i + mr 2 Mr 3 + rar 4 
X3 = M + m M + m ' ^ ' 

_ Mri + mr 2 + Mr 3 + mr 4 
x 4 = R = — . 5 

2M + 2m v 7 



We use the abbreviation x = {xi, x 4 } and r = {ri, r 4 }. 



A. The wave function 



The CG of the form 



1_ . , , 1 4 



G A (r) = exp{--f Ar} = exp{-- J2 A ij r i " r i) ( 6 ) 



is very popular in atomic and molecular physics [p]-|T0||. Here f stands for a one-row vector 
whose ith element is r$. The merit of this basis is that the matrix elements are analytically 
available and unlike other trial functions (for example, Hylleraas-type functions) one can 
relatively easily extend the basis for the case of more than three particles. The well-known 
defects of this basis are that it does not fulfill the cusp condition and its asymptotics does 
not follow the exponential falloff. This latter problem, especially for bound states, can be 
cured by taking linear combinations of adequately chosen CGs. 

The CG defined above is spherical and can thus describe systems with only L = orbital 
angular momentum. The usual way to account for the orbital motion in the case of L ^ 
is the partial-wave expansion. Because of the complexities arising from the evaluation of 
matrix elements this expansion gets very tedious for more than three particles. To avoid 
this difficulty the global vector representation [0] is used. In this approach, one defines a 
vector v as a linear combination of the relative coordinates: 
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v = ^w i r i , (7) 

i=l 

and the non-spherical part of the wave function is represented by a solid spherical harmonic 

yKLM(v)=V 2K+L Y LM (v)- (8) 
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The linear combination coefficients Ui are considered to be variational parameters and their 
optimal values are to be determined by the SVM as will be discussed later. The details and 
examples can be found in 

The calculation of the matrix elements for the space part of our basis function 

f KLM (u, A, r) = G A (r)y KLM (v) (9) 

is given in In the special case of K — the matrix elements can be written in much 
simpler form. This is shown in Appendix A. In the K ^ case, the CG is multiplied by 
a polynomial of the relative coordinates. In some cases this might be very useful, it can 
improve the short- distance behavior, for example, but this role can also be played by an 
appropriate superposition of the exponentials. We use K = in this paper. 

The translational invariance of the wave function is ensured by requiring that the param- 
eters A and u fulfill some special conditions. As is detailed in Appendix B, these conditions 
ensure that the motion of the center-of-mass is factorized in a product form. 

By combining the CG with the angular and spin parts, the full basis function takes the 
form 

§kLS = A{xsM s fKLM{u k , A k , r)}, (10) 

with an appropriate spin function Xsm S i where u k" is the index of the basis states and A is an 
antisymmetrizer for the identical fermions. In the positronium limit (er = 1) the Hamiltonian 
becomes invariant with respect to the interchange of positive and negative charges. Therefore 
the basis function should have a definite parity under the charge-permutation operator. See 
Appendix D for the details of the symmetry requirement on the wave function. For the 
special case with S — and Ms = in which two spins of positive charges and two electron 
spins are coupled to zero, respectively, the spin part of the wave function reads as 

xoo = i(| TTU) - I TUT) - I ITTI) + I UTT>)- (n) 

(Note that particles 1 and 3 are positive unit charges and particles 2 and 4 are electrons.) 
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Instead of optimizing the parameters of A it is more advantageous to rewrite Eq. (6) as 
ex p{ ~ r,) 2 -\H (12) 

i<j i 

The relationship between o;y,/3j and A is 

ay = -Aij (i ^ j), Pi = ^A ki , (13) 

k 

where acji (i < j) is assumed to be equal to Oy. There are two reasons to choose this form. 
The first is that in choosing in this way we deal with a correlation function between the 
particles i and j, while A^ has no such direct meaning and during the optimization it is 
more difficult to limit the numerical interval of A^ to be chosen from. Secondly, one can 
utilize this specific form to make the individual steps of the parameter selection very fast. 
By taking a look at the expressions of the matrix elements in Appendix A, it is clear that 
the main computational load is the calculation of the inverse and determinant of the matrix 
of the nonlinear parameters. The form in Eq. ( fP|) offers the possibility of the usage of the 
Sherman-Morrison formula to calculate these quantities, leading to a much faster function 
evaluation. The details of this step are given in Appendix C. 



B. Electric dipole transition rate 

In the positronium limit (a = 1) we deal with antiparticles and the electron-positron 
pair can annihilate. The lifetime of the first excited-state with L = 1 and negative parity 
is determined by both processes of annihilation and electric dipole transition to the ground 
state. The width r dipole for the latter process is calculated through the reduced transition 
probability B(E1) for the electric dipole operator = I]fc=i<?fc| r fc ~ R-|^i/i( r fc — R-) (/•* = 
-1, 0, 1) 

r dipo ie = ^(^) 3 5(^i;i--o + ), (14) 

with 
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B(El;l- ^0+) = J2\(00\D,\lM)\ 2 , (15) 
where E is the excitation energy of the first excited state. 

C. Annihilation rate 

The most dominant annihilation of the first excited-state of Ps2 is accompanied by the 
emission of two photons with energy of about 0.5 MeV each. The decay width T 27 for the 
annihilation can be estimated through the decay width of the para-positronium in spin- 
singlet state. This decay width has to be multiplied by the number N of positron-electron 
pairs which are in spin-singlet state in the Ps2- In the PS2 excited state we have four 
positron-electron pairs, among which the probability that the pair is in spin-singlet state is 
1/4 because the total spin of the first excited-state of Ps2 is zero, as will be shown later. 
(iVo = 4 x (1/4) = 1.) Therefore, to have an estimate for the decay due to the annihilation 



we can use the formula (2) of flO]] : 



r 27 = N T% (16) 



with 



r* = 4n(^) 2 hc(mn - r 2 )|*> = 4n(f) 4 hca^(6(r l2 )), (17) 

where the probability of finding an electron at the position of a positron, (5(r 12 )), is the 
expectation value of <5(xi) given in a.u., that is (5(ri 2 )) is equal to a^{^\5{ji — r 2 )| l I / ) with the 
Bohr radius ao- Roughly speaking, the lifetime is inversely proportional to the probability 
of finding an electron and a positron at the same position. 



D. The stochastic variational method 

To obtain very precise energy, one has to optimize the variational parameters u^i and 
A ki j of the trial function. The dimension of basis sets is typically between 100 and 1000, 
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and each basis state has nine nonlinear parameters. (See Appendix B.) The optimization of 
a function with a few thousands nonlinear parameters cannot be done efficiently by using a 
deterministic optimization method, since this could entail the complete reconstruction of the 
Hamiltonian matrix and diagonalization every time when some of the nonlinear parameters 
are altered. Moreover, the deterministic search for the optimal value of such a large number 
of parameters is likely to get trapped in a local minimum. 

A procedure based on the stochastic search for the best set of nonlinear parameters can 
be programmed efficiently 0,|lTJ and is capable of achieving highly accurate results for most 



few-body systems p|,^j7|JT^| . The essence of the strategy can be summarized as follows: Let 
{ui, Ai] be the nonlinear parameters of the ith basis function out of the set of K such basis 
functions. Then the procedure is 

(1) A succession of different sets of ({it,], A]}, {w™ s ,v4™ s }) are generated randomly. 

(2) By solving the eigenvalue problem, the corresponding energies (Ef, E™ 3 ) are deter- 
mined. 

(3) The parameter set {u™, A™} which produces the lowest energy is then used to replace 
the existing {ui, Ai} set. 

(4) The procedure cycles through the different parameter sets ({ui,Ai},i = 1,...,K), suc- 
cessively choosing different sets to minimize the energy until convergence is reached. 

The essential reason motivating this strategy is the need to sample different sets of nonlinear 
parameters as fast as possible. The main advantage is that it is not necessary to recompute 
the complete Hamiltonian nor it is necessary to solve the generalized eigenvalue problem 
from scratch each time a new parameter set is generated. By changing the elements of 
parameter set for each basis function individually, it is necessary to recompute only one row 
(column) of the Hamiltonian and overlap matrices each time the parameter set {ui, Ai} is 
changed. Furthermore, the solution of the generalized eigenvalue problem is also expedited 
since the Hamiltonian matrix is already diagonal apart from one row and one column. 
A similar strategy to the above was used when adding additional terms to the basis. 
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The speed of the calculation can be further increased if one changes the nonlinear pa- 
rameters Ai in a special way. This is described in Appendix C. 

The above way of finding the best parameters is certainly very restricted. Even this 
simple method gives very accurate energies. More sophisticated technique may give better 
results in a smaller basis size. 

III. RESULTS 

The results of calculations for the ground state of HPs and Ps 2 and the first excited-state 
of Ps 2 are reported in this section. The ground states of HPs and Ps 2 have already been 
subject to intensive calculations and some of the results obtained before for these systems 
serve as validation of the SVM. The calculation of the properties of the excited state of the 
Ps 2 is the primary focus of this paper. We have previously reported the energy of the ground 
state of the Ps 2 and predicted the existence of an excited state of this molecule. This paper 
reports considerably improved energies by further optimization of the nonlinear parameters 
of the basis. The further optimization and the increase of the basis dimension has produced 
an improved wave function and we present different properties of these systems by using that 
wave function. We show the convergence of the binding energies and various expectation 
values as a function of the dimension of the basis. The results in the tables are shown for 
the basis dimensions of K = 100,200,400,800,1200,1600. The basis has been subject of 
intensive optimizations at these dimensions. Once the optimization at a given basis size has 
been finished, new basis states have been added (each of them has been selected amongst 
hundreds of random candidates) to reach the next basis size where the optimization is started 
again. While the pattern of convergence is a very useful information about the accuracy of 
the results, one has to keep in mind that this can be distorted by many extraneous factors. 
This is because one cannot guarantee that the quality of these optimizations is the same. 
We expect that the stochastic selection of the basis is close to be the optimal choice for 
lower dimensions, but for large dimensions (K = 1200, 1600) the procedure becomes more 
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time consuming and we have less chance to find the optimal parameters. 



A. Hydrogen Positride, HPs 

The boundness of the exotic molecule, HPs, has been known theoretically for many 



years [13| and it has recently been created and observed in collisions between positrons and 



methane [II]. The investigation of the stability of positronic atoms has been attracting much 
attention because positrons can be used as a tool for positron-annihilation spectroscopy in 
condensed matter physics. The HPs molecule is the simplest but ideal hydride to test the 
SVM. It is also very intriguing to see the difference between the properties of Ps2 and HPs. 

The energy calculated by SVM and by other methods are shown in Table I. The proton 
mass is taken to be infinite. The two electrons are assumed to be in spin-singlet state. The 
spin states of the proton and the positron can be taken arbitrary. Our result, already at 
the dimension of K = 200, is better than the previous calculations. The increase of the 
basis size improves the energy further. The need of improved accuracy can be clearly seen 
in Table II, where various expectation values are listed. The expectation value (rg_ e _), for 
example, is much less accurate than the energy and it is considerably improved beyond the 
dimension K = 200. 

One can compare the expectation values of the separation distances of the particles in 
the HPs to those in the H and Ps atoms. The average electron-positron distance is 3.48 ao 
in HPs, which is slightly different from that in the positronium atom (3 a ). The average 
electron-proton distance in HPs and H is considerably different (2.31 a and 1.5 a ). The 
average distance between the two positive charges (3.66 ao) is much larger than that in the 
H2 molecule (1.41 ao). 

The correlation function defined by 

C{r) = {*\5(r i -r j -r)\*) (18) 

gives more detailed information on a system than just various average distances. This 
quantity can be calculated by using Eqs. (34) and (35). For the spherical wave function 
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with L = 0, C(r) is a function of r, that is, the monopole density, and for the L = 1 wave 
function, it consists of two terms of monopole and quadrupole densities. Figure 1 displays 
r 2 C(r) for various pairs of the constituents of HPs. The two electrons are attracted by 
the proton, but the proton-electron correlation function is much broader than that in the 
H atom, while they are separated with its maximum density being at about 2.8 a.u. The 
positron moves furthest from the proton and has a peak density around at 2.6 a.u. from the 
electron. 

The 27 annihilation rate, calculated from Eq. (16) with Nq = 2 x (1/4) = 1/2 and 
(5 e + e -) of Table II, is found to be Y 2l = 2.4722 ns -1 , improving the previous estimates |21 
by about 0.5%. 



B. Positronium molecule, PS2: Ground state 

The energies by SVM are compared to the best previous results in Table III. The result of 
SVM, again, already at the dimension of K = 200, is better than the energy of the previous 
calculations. The increase of the basis size improves the accuracy and the virial factor 
|1 + {V)/(2(T))\ becomes as small as 0.3 x 10~ 9 , improving the previously best calculation 
by more than 4 order of magnitude. 

The average electron-positron distance is 4.487 ao, which is about 1.5 times larger than 
in the positronium atom. The 27 annihilation rate calculated from Eq. (16) by using (5 e + e -) 
of Table II is found to be T2 7 = 4.470 ns -1 . 

The electron-electron and the electron-positron correlation functions are compared in 
Fig. 2. The peak position of the electron-electron correlation function is shifted to larger 
distances compared to the one of the electron-positron correlation function. The electron- 
positron correlation function in PS2 has much broader distribution than the corresponding 
function in a Ps atom. 
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C. Positronium molecule, PS2: First excited state 



In our previous paper we have predicted the existence of the first excited-state of the 
Ps 2 molecule. This is a unique bound-state which cannot decay into two Ps atoms due to 
the Pauli principle. The spin of this state is S = and the orbital angular momentum is 
L = 1 with negative parity. In this spin state, the Ps 2 molecule can dissociate into two Ps 
atoms (bosons) only if the relative orbital angular momentum is even. Consequently, the 
Ps2 molecule with L = 1 and negative parity cannot decay into the ground states of two 
Ps atoms (Ps(L = 0)+Ps(L = 0)). The energy of this Ps 2 (L = 1) state (E = -0.334408 
a.u., see Table III) is lower than the energy of the relevant threshold (—0.3125 a.u.), and 
this state is therefore stable against the autodissociation into Ps(L = 0)+Ps(L = 1). The 
binding energy of this state is 0.5961 eV, which is by about 40% more than that of the 
ground state of Ps 2 (0.4355 eV). 

We have shown in that the bound excited state is essentially a system where two Ps 
atoms, one in its ground state and the other in its first excited P state, are weakly coupled. 
The expectation value of the average electron-positron distance shown in Table V supports 
this picture: The value of 7.57 ao in the excited state is 15 % larger than the average (6.5 
ao) of the electron-positron distances in the L = ground state of the Ps atom (3 ao) and 
the L = 1 excited state of the Ps atom (10 a ). We can also estimate the root-mean-square 
distance d = J (x 3 • x 3 ) between the two atoms by 

j2 = ^ ^ + r, _ £3+^ = 1 ^ + ^ _ ^ ^ (u) 

The symmetry properties of the Ps 2 wave function are used to obtain the second equality. 
Using the values of Tables IV and V yields d = 6.93 a.u. for the L = 1 excited state and 
d = 4.82 a.u. for the L = ground state. 

Figure 3 displays the electron-electron and electron-positron correlation functions. As 
mentioned before, the correlation function for the L = 1 state consists of the monopole and 
quadrupole densities and their shapes depend on the magnetic quantum number M of the 



12 



wave function. Of course the M-dependence of the shapes is not independent of each other 
but is determined by the Clebsch-Gordan coefficient. See Eq. (34). The quadrupole density- 
is contributed only from the P-wave for the electron-positron relative motion, while the 
monopole density is contributed by both S- and P- waves. Figure 3(a) plots the correlation 
functions for M = and Fig. 3(b) the correlation functions for M — 1. As the correlation 
function is axially symmetric around the z axis and has a reflection symmetry with respect 
to the xy plane, the correlation function sliced on the xz plane is drawn as a function 
of x (x > 0), z (z > 0). The electron-electron correlation function has a peak at the point 
corresponding to the average distance of 7.57 a.u. The electron-positron correlation function 
has two peaks reflecting the fact that the basic structure of the second bound-state is a 
weakly coupled system of a Ps atom in the L = state and another Ps atom in the L — 1 
spatially extended state. The peak located at a larger distance from the origin is due to the 
P-wave component of the Ps2 molecule. 

By using the obtained value for (5(r 12 )) in Eq. (16), the lifetime due to the annihilation 
is estimated to be 0.44 ns. This is about twice of the lifetime of the Ps 2 ground state. The 
B(El) value is calculated to be B(E1) = 0.87e 2 aQ. By combining this value with the dipole 
transition energy of 4.94 eV, the lifetime due to the electric dipole transition has been found 
to be 2.1 ns. The branching of the electric dipole transition is thus about 17 % of the total 
decay rate. Therefore, both branches contribute to the decay of the excited state of the Ps 2 
molecule. Its lifetime is finally estimated to be about 0.37 ns. The excitation energy of 4.94 
eV found for the Ps 2 is different by 0.16 eV from the corresponding excitation energy (5.10 
eV) of a Ps atom. This difference seems to be large enough to detect its existence, e.g. in 
the photon absorption spectrum of the positronium gas. 

IV. SUMMARY 

We have used the Correlated Gaussians combined with the angular functions which are 
specified by the global vector. Nonlinear parameters of the bases have been determined 
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by the stochastic variational method. We have considerably improved the results of the 
previous calculations for the ground state of HPs and Ps2- In addition, we have calculated 
various expectation values, correlation functions and other properties of the excited state of 
the Ps2 molecule. 

The excited state of the PS2 molecule has the orbital angular momentum L = 1, the 
spin S = 0, and negative parity. The excitation energy of the state is 4.941 eV, 0.596 eV 
below the threshold of Ps(L = 0)+Ps(L = 1). Though this state is in the continuum of 
Ps(L = 0)+Ps(L = 0) channel, it is stable against autodissociation into that channel because 
of the Pauli principle. The main decay mode of this state is the annihilation emitting two 
photons of about 0.5 MeV each which, except for the tiny binding-energy difference, is equal 
to the photon energies of Ps atoms. The annihilation decay mode is not useful to confirm 
experimentally the existence of the Ps 2 molecule. We have discussed a unique decay mode 
of the excited state, the electric dipole transition to the ground state. The lifetime due to 
the electric dipole transition has been calculated to be 2.1 ns, while the lifetime due to the 
annihilation is 0.44 ns. The electric dipole transition can be used as a signal for experimental 
confirmation of the Ps2 molecule. 

This work was supported by OTKA grant No. T17298 (Hungary) and Grant-in-Aid 
for International Scientific Research (Joint Research) (No. 08044065) of the Ministry of 
Education, Science and Culture (Japan). The authors are grateful for the use of RIKEN's 
computer facility which made possible most of the calculations. 

APPENDIX A: EVALUATION OF MATRIX ELEMENTS 

In this appendix the matrix elements of the spatial part of the basis functions are given. 
The method of calculation of these analytical expressions is detailed in Refs. [0,^2] • The 
main aim of this section is to convince the reader that the formulae are particularly simple 
for the case of K = 0. The extension to a general iV-body system is straightforward so that 
we assume that the system contains N particles. 
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The basic idea of the calculation of the matrix elements is the usage of the generating 
function g: 

g(s;A, r) = exp ( - ^fAr + sr). (20) 
In the special case of K = 0, Eq. (|J) is obtained from g by 
foLMiu, A, r) = e-^ Ar y 0L Ai(v) 



^ J Y LM (e) ( J^<?(Aeu; A, r) j de, (21) 



A=0,e=l 



with 



To abbreviate the expression for the matrix elements we introduce the following notation 

(f'\0\f) = (f 0LM (u',A', r )\O\f 0LM (u,A,v)), (23) 

where O stands for the unity, kinetic or potential energy operators. The operators considered 
here are rotational invariant and thus the matrix elements are diagonal in LM. Note the 
prime on / which is a reminder that the parameters in the ket and the bra may be different. 

The use of Eq. fl2~ID in Eq. ( p3|) leads to an expression that the matrix element is derived 
from that between the generating functions, which becomes a function of parameters A, e, A' 
and e'. Here the matrix element between the generating functions can be obtained easily by 
using the expression 



1 detA ' 



and its extended formulae. After a power series expansion of the matrix element between the 
generating functions in terms of A, e, A' and e', the derivative and the integration prescribed 
in Eq. ( PT| ) can be carried out straightforwardly |7|,|22| . 
The overlap of the trial functions is given by 

^ = (^)' B ^- (25) 

15 



The kinetic energy is expressed by 



(f'\T - T cm \f) = y( R + L ®P~ l ) (/'!/>■ (26) 



The matrix elements of a central potential reads as 

L\ /77^ 

l{C,n)—r 

where the integral over the radial form of the potential is expressed with use of Hermite 
polynomials 



(f\V(\*i ~ r;l)l/> = (f'\f) E ^rA^ (t^)" (27) 



1 /"°° /2 

J (^) = ^ { 2n + iy. l V^-*y x2 Hi(x)H 2n+1 (x)dx. (28) 
The definitions of the constants in the above expressions are 



B = A + A', p = u'B~ l u, p = p - -77'. 



R = STriB^A'AA), Q = 2u 'B' 1 AAA'B^u. 

c- 1 = vfitiB^w^, 7 = cw^3)B~ l u, 7' = c«^)S _1 u', (29) 

where the N x N symmetric matrix A is defined by T — T cm = (1/2) J2i,j ^ijPi • Pj an d 
is an N x 1 one-column matrix defined by 

w { k ij) = 5 ki -5 kj (k = l,...,N). (30) 

The integral in Eq. (28) can be analytically evaluated for several potentials, including 
Coulomb, exponential or Gaussian potentials. The numerical evaluation for a general po- 
tential is a simple matter and one tabulates J(c, n) for the necessary values of c. For power 
law potentials V(r) = r fc , for example, including the Coulomb interaction, the c-dependence 
of the integral I{c,n) is factored out: 

/2 x fc/2 

I(c,n) = (-j I n (k), (31) 
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where the remaining integral can be carried out and expressible in terms of the Gamma 
function: 

k + 3^ 



i n ( i \mo2n- 2m+l 

i (k) = — y { ' 

^ ^ m\(2n-2m + l) 



m + 



In particular, for the Coulomb force (k = —1) we get 

^4 (-l) n 



4-1 



7r (2n + l)n! 

The correlation function is calculated through the equation 



where 



(f|5(r J -r,-r)|/) = (f|/)cle-^ 2// '' 



n=0 



x ^C LnK (LM2 K 0|LM)F 2K0 (r) 

K=0 



C 



LriK 



(-1) k (2k - (2L - 2k)!(2L + 2k + 1)! /4k + 1 
tt(2L - 1)!!2 L+1 /2(L - n)!K!(n - k)!(2w + 2k + 1)!! V 2L + 1 ' 



(32) 



(33) 



(34) 



(35) 



APPENDIX B: SEPARATION OF CENTER-OF-MASS MOTION 



The transformation between the relative and single-particle coordinates, given by Eqs. 
-(I), can be defined by the matrix: 



U 



1 


-1 














1 


-1 


M 


m 


— M 


— m 


m+M 


m+M 


m+M 


m+M 


M 


m 


M 


m 



\ 



( 



u 



-1 



\ 2m+2M 2m+2M 2m+2M 2m+2M J 



m+M 
-M 



o I 1 



1 1 

m+M u 2 x 





V 



m+M 2 
m+M ~2 1 / 



(36) 



The transformation between the relative and the single-particle coordinates is given by 

x = Ur, r = U _1 x. (37) 

Here r and x are column vectors containing (ri, ...,r4) and (xi, ...,X4). By this transforma- 
tion one can express the CG of the single-particle coordinates by the relative coordinates: 
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G a (t) = exp{--fAr} = exp{-^xAx} = G A (x), A = U- 1 AU~ 1 . (38) 

The parameters Ajvi = A^, (i = 1,...,N — 1) connect the relative and center-of-mass 
variables, and give rise to an undesirable center-of-mass dependence of the wave function. 
To have a translational invariant basis, we require that 

Am = 0, A NN = c, 

that is, 

N N N N 

EE^ 1 = (< = i,...,iv-i), EE^ = c . ( 39 ) 

3=1 k=l 3=1 k=l 

where c is an arbitrary, positive constant common for each basis function. The second 
condition assures the finite norm of the basis function. By this requirement the relative and 
center-of-mass motion is separated in the exponential part of the basis function. 

To remove the center-of-mass contamination from the angular part, let us express the 
global vectors v in terms of relative coordinates: 

N N N 

V = J2 U i r i = J2 U iJ2 U ik*k- ( 4 °) 
i=l i=l k=l 

This identity shows that by requiring 

N N 

E«i c/ w=E«i = °. ( 41 ) 
i=i i=i 

the global vector becomes translationally invariant. 

By fulfilling Eqs. (39) and (41) the basis is free from any problems with the center-of- 
mass motion. These conditions fix iV + 1 nonlinear parameters among iV(iV + l)/2 + iV = 
N(N+3)/2 parameters. For N = 4 there remain nine free parameters for each basis function. 



APPENDIX C: SHERMAN-MORRISON FORMULA 



As it is shown in Appendix A, the calculation of the matrix elements requires the evalu- 
ation of the determinant and inverse of the matrix B. In the SVM process we probe many 
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random trials with different matrices. Let us assume that we change the matrix A of non- 
linear parameters in such a way that we change the parameter (i ^ j) of the relative 
motion between particles % and j to + A but keep all other matrix elements unchanged. 
This is certainly a very restricted way, but in this case the computer time required for the 
evaluation of the matrix elements tremendously decreases. This change of produces the 
following changes in the matrix A (see Eq. (13)): 

Aij > Aij A, Aji > Aji A, An > An + A, Ajj > Ajj + A. (42) 

It is easy to see that this change does not violates the conditions of Eq. (39). Thus the 
wave function with this modification is still translational invariant. The above change in the 
matrix A can be simply expressed by using the vector defined in Eq. (30) as follows: 

A ^ A - Xw {lj) w&\ (43) 

Note that w^w^ is an N x N matrix, whereas w^w^ is just a number. As B is equal 
to A + A', the above change leads to the following modification of B, 

B -> B- \w {ij) u0l (44) 

To calculate the inverse and determinant of the above special form, the Sherman- 
Morrison formula can be used: 

(B- Xw^w&A" 1 = B- 1 + — ^ B- l w {lj) w^i)B-\ (45) 

v ' 1 - Xw^B^w^) 

and 

det (B - Xw iij) u0^ = (l - X^B^w^) detS. (46) 

The advantage of this formulae is apparent: By knowing B^ 1 and det-B one can easily 
calculate the right-hand sides of the equations, and the A dependence is given in a very 
simple form. For example, w^B^w^ simply reduces to (B~ l )u + (B~ l )jj — 2(B~ 1 ) ij . 
Likewise, B^ 1 w^w ( ~ i ^B~ 1 can also be easily evaluated. To change A, therefore there is no 
need for the evaluation of inverses and determinants (which would require N 3 operations) 
but we get the desired results by a simple multiplication and division. 
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APPENDIX D: SYMMETRIZATION OF WAVE FUNCTIONS 



Antisymmetrization 



The antisymmetrizer A is defined as 

A=- 



■y n P 



in 



(47) 



V i=l 



where the operator Vi changes the indices of identical particles according to the permutation 
■■■p t N ) of the numbers (1, 2, N), and is the phase of the permutation. The effect of 
this operator on the set of the position vectors (ri, rjv) is 



V i (r 1 ,...,r N ) = (r pl ,...,r pV ) 



(48) 



By representing the permutations by the matrix 



( C i)kj = 1 if i = Pk and ( C i)ki = otherwise, 



kj 



(49) 



(for example, the permutation (3 12 4) is represented by 



^0 10^ 



C 



10 
10 
1 



(50) 



while for (1 2 3 4) C is a unit matrix), the effect of the permutation operator on the single- 
particle coordinates reads as 



ViV = Qr. 



(51) 



By using Eqs. (37) and (51) the permutation of the relative coordinates is expressible as 



P.x = fix with P^UCJJ- 1 



(52) 



The CGs, after permutation, take the form: 
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P,G A (r) = G diACi (r) = Gp iAPi (x) = G 



x . 



(53) 



In the spin space the permutation operator interchanges the indices of the single-particle 
spin functions and can be easily evaluated. As a result, the matrix element of any spin- 
independent operator O which is symmetrical with respect to the permutation of identical 
particle coordinates can be written in the following form: 



(A{xsM s fKLM(u', A', r)}\0\A{xsM s fKLM(u, A, r)}) 

n p 

£ Ci{fKLM(u', A', v)\0\f KLM (C t u, C t AC u r)>, 



i=l 



where the coefficients c,- have the form 



Ci — £i\XSM s \Pi\XSMs) ■ 



(54) 



(55) 



Since the antisymmetrizer is a projector onto an antisymmetric state, only ket (or bra) 
function needs to be antisymmetrized. 

The particular value of the coefficient q depends only on the spin function of the system. 
In the case of Ps2 two positrons must be antisymmetrized and likewise two electrons must 
be in antisymmetric states. Therefore, the antisymmetrizer for this system is given by 
A — (1 — -P13XI — P24), where Pij is the transposition of particle labels i and j. Thus A has 
four permutations (n p = 4) and we can identify 

Pi = l, V 2 = Pi 3 , V 3 = P 2i , P 4 = P 13 P 24 . (56) 
The corresponding phases are E\ = 1, e 2 = — 1, £3 = — 1, £4 = 1 and the matrices C are given 



as follows 



10 
10 
10 
1 



c 2 



^0 10^ 



10 
10 
1 



^ 1 ^ 



1 
10 
10 



Ca 



^0 10^ 



1 
10 
10 



(57) 



The spin function xoo of Eq. QTTj) is antisymmetric in both of the positron spin coor- 
dinates and the electron spin coordinates. Thus the spin matrix element (xsM s \Pi\XsM s ) 
turns out to be equal to £j and we have c\ — c 2 = C3 = C4 = 1. 
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Charge symmetry 



The Hamiltonian H for Ps 2 has charge-exchange symmetry, that is, it is invariant under 
the exchange of the positive and negative charges: Letting P denote the charge-permutation 
operator, we have 

HPip = PHi/j = EPtp. (58) 

Therefore, the non-degenerate eigenstate of the Hamiltonian is also the eigenstate of the 
charge-permutation operator. In the Ps2 the ground state is even (tt = +1) under P, while 
the L = 1 excited state turns out to be odd (n = —1). 

Consider the case of Ps (e + e~). This system is represented by the coordinate (r 2 — r^). 
The charge permutation is thus equivalent to the parity operation. Since the parity is 
(— 1) L for the state with orbital angular momentum L, the eigenvalue of charge-permutation 
operator is also (— 1) L . The signs of xi and x 2 change with respect to the charge permutation 
-P12-P34, while x 3 does not. Assume that the Ps 2 has partial waves l±, / 2 and / 3 corresponding 
to the motion described with x x , x 2 and x 3 , respectively. When charges are permutated, 
the wave function ip become (— l)' 1 (— l) l2 ip. Then Pip = ip for the S state with li — l 2 — 0, 
while Pip = —ip for the P state with li = and Z 2 = 1. 

The non-vanishing matrix element of the electric dipole transition supports that the first 
excited P-state is odd under the charge permutation. This is because the electric dipole 
operator D has the following form except for the constant: 

D = e(n - R) - e(r 2 - R) + e(r 3 - R) - e(r 4 - R), (59) 

which changes sign under the charge permutation. Therefore, if the excited P-state is even 
under the charge permutation, then the electric dipole matrix element between the P state 
and the ground state would identically vanish. 

The charge-permutation operator P is given by Pi 2 P 3 4 or Pi4P 32 . When the wave function 
■if) is already antisymmetrized for two positrons and for two electrons, then we can see that 
both operators give the same effect. To understand this we use the following identity 
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PuP&*l> = {Pl2PM) 2 PlAP?,2{PlzP2±) 2 i> = P^PsP^ = ^12^34^- 



(60) 



Thus the basis function for the PS2 molecule with definite charge-permutation symmetry is 
given by operating with the following operator C on the function: 

1 



C = -^(1 + nP 12 P u )(l - P 13 )(l - P 24 ) 



(61) 



The evaluation of matrix elements between the states with odd charge symmetry can be 



done in a similar manner to the previous subsection by extending Eqs. (p4|) and (|55). The 
antisymmetrizer A is now replaced with C = (l/v8) Yn=i e iPu where new Vi are defined by 



P* — p2-p4; ^6 — p2-p4p3) Pi — PyiPmPiA, Pr — P 2pU -Pl3p 



12- r 34- r 13- r 24; 



(62) 



and the corresponding phases are e$ = —1, eq = 1, £7 = 1, £s = — 1- The matrices Cj 
corresponding to Pj are given below: 



^0100^ 



1000 
0001 
0010 



c fi 



^0001^ 



1000 
0100 
0010 



^0100^ 



0010 
0001 
1000 



^0001^ 



0010 
0100 
1000 



• (63) 



It is easy to evaluate the coefficients q. For the spin function xoo of Eq. ([□])> we get 

c 5 = c 6 = c 7 = c 8 = ~~ 1- 
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TABLES 

TABLE I. Comparision of the results of different calculations for the ground-state energy of 
HPs. The energy is given in atomic units. 



Method 



Reference 



Energy 



SVM (K = 100) 
SVM {K = 200) 
SVM (K = 400) 
SVM {K = 800) 
SVM (K = 1200) 
SVM (K = 1600) 

Hylleraas configuration interaction 

Exponential trial functions 

Diffusion Monte Carlo 

Diffusion Monte Carlo 

Correlated Gaussian basis {K = 200) 



Present 
Present 
Present 
Present 
Present 
Present 



15 



|17| 



[18| 



E3 



-0.7891013600 

-0.7891810473 

-0.7891924458 

-0.7891958706 

-0.7891964226 

-0.7891965536 

-0.7842 

-0.7889 

-0.7891 ±0.002 
-0.789175 ± 0.00001 
-0.7891794 
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TABLE II. Expectation values of various quantities for HPs. Atomic units are used. K is the 
basis dimension. 









E 


-(V)/(2(T)) 






K 




100 


-0 7891013600 


1 00001 






K 

J. V 




900 


—0 7891 81 0473 


1 000003 






K 




400 


—0 7891 924458 


1 000002 






K 




800 


-0.7891958706 


1.0000007 






K 




1200 


-0.7891964226 


1.0000004 






K 




1600 


-0.7891965536 


1.0000003 












(r 4 ) 


V e+e-l 


(r 4 ) 

V e-'pl 


V e+p/ 


K 




100 


515.42669 


525.13203 


193.45055 


504.56556 


K 




900 


594 98363 


5^1 94495 


1 97 60909 


51 3 48089 


K 




400 


527 33506 


532 59188 


1 98 63996 


515 59169 


K 




800 


527.88970 


532.94707 


198.88278 


516.06972 


K 




1200 


527.94660 


532.98328 


198.90610 


516.11702 


K 




1600 


527.96159 


532.99639 


198.91176 


516.13646 










V e+e-/ 


(r 3 ) 


\ e+p/ 


K 




100 


83.599992 


83.792382 


34.789685 


84.226327 


K 




900 


84 S37498 


84 949983 


35 1 90409 


84 91 1 659 


K 




400 


84 507962 


84 347282 


35 195647 


85 0641 1 2 


K 




800 


84.544707 


84.369687 


35.211858 


85.094386 


K 




1200 


84.548681 


84.372106 


35.213444 


85.097517 


K 




1600 


84.549852 


84.372949 


35.213895 


85.098746 








/ r 2 v 
\ e e I 


(ri.e-) 


/ r 2 v 
\ e pi 


(r 2 e+P ) 


K 




100 


15.803193 


15.542251 


7.7797451 


16.188998 


K 




200 


15.860043 


15.575673 


7.8062352 


16.241186 


K 




400 


15.872464 


15.582575 


7.8117324 


16.252128 
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K 




800 


15.874993 


15.584009 


7.8128668 


16.254178 


K 




1200 


15.875286 


15.584176 


7.8129800 


16.254399 


K 




1600 


15.875377 


15.584230 


7.8130152 


16.254480 








( r e-e-) 


( r e+e-) 


( r e-p) 


( r e+p) 


K 




100 


3.5700072 


3.4777333 


2.3092381 


3.6573544 


K 




200 


3.5738023 


3.4797561 


2.3110943 


3.6607696 


K 




400 


3.5745993 


3.4801765 


2.3114423 


3.6614669 


K 




800 


3.5747568 


3.4802575 


2.3115152 


3.6616016 


K 




1200 


3.5747763 


3.4802676 


2.3115221 


3.6616167 


K 




1600 


3.5747825 


3.4802707 


2.3115245 


3.6616220 








(C-V) 

' e e ' 




N e p 1 


KL) 

N e+p' 


K 




100 


0.37072021 


0.41851818 


0.72973620 


0.34760250 


K 




200 


0.37058889 


0.41850815 


0.72971467 


0.34749891 


K 




400 


0.37056069 


0.41849668 


0.72970918 


0.34746907 


K 




800 


0.37055594 


0.41849614 


0.72970858 


0.34746293 


K 




1200 


0.37055519 


0.41849601 


0.72970874 


0.34746209 


K 




1600 


0.37055494 


0.41849596 


0.72970869 


0.34746180 








<v-V> 

N e e ' 




<v- 2 „> 

N e p 1 


N e^pi 


K 




100 


0.21426165 


0.34877458 


1.2059515 


0.17234727 


K 




200 


0.21396622 


0.34911573 


1.2069510 


0.17221620 


K 




400 


0.21392019 


0.34912443 


1.2070112 


0.17217310 


K 




800 


0.21391300 


0.34914011 


1.2070561 


0.17216589 


K 




1200 


0.21391137 


0.34914210 


1.2070629 


0.17216413 


K 




1600 


0.21391064 


0.34914275 


1.2070632 


0.17216372 








^4 e~ ' T eae+ ) 


( T e+e- - r e+e-) 


(r --r -) 


(^6+ 'Tpe~ ) 


K 




100 


7.9015967 


7.6406546 


-0.12185159 


4.2132458 


K 




200 


7.9300217 


7.6456510 


-0.12378653 


4.2358745 
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K 




400 


7.9362320 


7.6463425 


-0.12449962 


4.2406428 


K 




800 


7.9374963 


7.6465132 


-0.12462952 


4.2415176 


K 




1200 


7.9376432 


7.6465328 


-0.12466320 


4.2416014 


K 




1600 


7.9376883 


7.6465421 


-0.12467313 


4.2416325 










"< V e + > 




(V e + • V e -) 


K 




100 


0.65224870 


0.27367198 


-0.043864431 


0.11701815 


K 




200 


0.65232846 


0.27369666 


-0.043999455 


0.11707408 


K 




400 


0.65234077 


0.27369750 


-0.044052593 


0.11707637 


K 




800 


0.65234481 


0.27369980 


-0.044060768 


0.11707718 


K 




1200 


0.652345728 


0.27370016 


-0.044063957 


0.11707760 


K 




1600 


0.652345903 


0.27370022 


-0.044064366 


0.11707739 








(<W> 


(<5 e + e -) 


(S e - P ) 


(^e+p) 


K 




100 


0.0047127 


0.0236658 


0.1717649 


0.0017964 


K 




200 


0.0047873 


0.0242912 


0.1758767 


0.0016985 


K 




400 


0.0044178 


0.0243887 


0.1761969 


0.0016542 


K 




800 


0.0043895 


0.0244224 


0.1768711 


0.0016440 


K 




1200 


0.0043889 


0.0244583 


0.1771854 


0.0016386 


K 




1600 


0.0043867 


0.0244611 


0.1771862 


0.00163857 
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TABLE III. Total energies of the Ps2 molecule in atomic units. K is the basis dimension. 



Method 


Ps 2 (L = 0) 


Ps 2 (L = 1) 


SVM (K = 100) 


-0.516000069 


-0.334376975 


SVM (K = 200) 


-0.516003119 


-0.334405047 


SVM (K = 400) 


-0.516003666 


-0.334407561 


SVM (K = 800) 


-0.516003778 


-0.334408177 


SVM (K = 1200) 


-0.5160037869 


-0.334408234 


SVM (K = 1600) 


-0.516003789058 


-0.3344082658 


Ref. Q (K = 200) 


-0.5160024 




QMC @ 


-0.51601±0.00001 
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TABLE IV. Expectation values of various quantities for the ground state of Ps2- Atomic 



units are used. The positrons are labeled 1 and 3 and the electrons are 2 and 4. Because of the 
charge-permutation symmetry, e.g. (r\2) = (ru) = (r^z) = ( r 34)- K is the basis dimension. 
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vw 


/V \ 
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K = AOO -0.25800178 O.lxlO" 6 

K = 800 -0.25800188 0.2xl0~ 7 

K = 1200 —0.25800188 0.4xl0~ 8 

if = 1600 -0.258001894 0.3xl0~ 9 



TABLE V. Expectation values of various quantities for the excited state of Ps2- Atomic units 
are used. See the caption of Table IV. 
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K = 800 -0.1672038 0.8xl0~ 6 
K = 1200 -0.1672039 0.5xl0~ 6 
K = 1600 -0.16720401 0.36xl0~ 6 
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FIGURES 
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Fig. 1 : The correlation functions r 2 C (r) for various pairs of the constituents of the hydrogen 
positride HPs. For the sake of comparison, the electron-proton correlation function of 
the H atom is also drawn. 
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2 : The correlation functions r 2 C(v) for the ground state of the Ps 2 molecule. The solid 
curve denotes the electron-electron correlation and the dashed curve the electron- 
positron correlation. For the sake of comparison, the electron-positron correlation 
function for the Ps atom is drawn by the dotted curve. 
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Fig. 3 : The correlation functions r 2 C(r) (r = (x, 0, z)), multiplied by one thousand, for 
the bound excited-state of the PS2 molecule which has the orbital angular momentum 
L = 1, the spin 5 = 0, and negative parity. The magnetic quantum number M is set 
equal to for (a) and to 1 for (b). Plotted on the xz plane are the contour maps of 
the correlation function. 
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